Quantum phases in an optical lattice 
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We present the zero-temperature phase diagram of bosonic atoms in an optical lattice, using 
two diflFerent mean-field approaches. The phase diagram consists of various insulating phases and a 
superfluid phase. We explore the nature of the insulating phase by calculating both the quasiparticle 
and quasihole dispersion relation. We also determine the parameters of our single band Bose- 
Hubbard model in terms of the microscopic parameters of the atoms in the optical lattice. 
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I. INTRODUCTION 

Using the interference pattern of intersecting laser 
beams one can create a periodic potential for atoms, 
which is known as an optical lattice [Qj2j. Because one 
can confine atoms at separate lattice sites, one can accu- 
rately control the interaction between the atoms. This 
makes the optical lattice an important tool in spec- 
troscopy, laser cooling and quantum computing Q . In 
the following we study some of the many-body aspects of 
such a lattice and in particular Bose-Einstein condensa- 
tion of atoms in an optical lattice. In contrast with the 
existing Bose-Einstein condensation experiments in an 
harmonic trap, the quantum depletion of the condensate 
in the case of Bose-Einstein condenstation in a lattice 
can be very large. We can therefore expect interesting 
features. 

If we assume that the atoms are cooled to within the 
lowest Bloch band of the periodic potential, Jaksch et 
al. have shown that we can describe the behavior of the 
atoms in an optical lattice with the Bose-Hubbard hamil- 
tonian 



H 



(1) 



where the sum in the first term on the right-hand side 
is restricted to nearest neighbours and and are the 
creation and annihilation operators of an atom at site i 
respectively. The parameter t is the hopping parameter 
and U is the interaction strength, which we always as- 
sume to be positive in the following. The term involving 
the chemical potential /i is added because we perform our 
calculations in the grand-canonical ensemble. 

Qualitatively we expect that when there is an integer 
number of particles at each site i and t U, the in- 
teraction between the particles will make it energetically 
unfavourable for a particle to move from one site to an- 
other. In this situation the gas is in what is known as the 
Mott insulator phase Q| . However, if we add in this phase 
a particle to the system, this particle will only receive a 



small energetic penalty when it moves, because its inter- 
action energy is the same on each site. For this reason, 
a gas with a non-integer number of bosons at each site 
will be in a superfluid phase at zero-temperature. This 
expectation has been shown to be correct using Quan- 
tum Monte Carlo calculations and several mean-field 
approaches [^|j8|J^. In particular, Ref. Q numerically de- 
termines interesting features of cold bosonic atoms in an 
inhomogeneous optical lattice. In this paper, we give a 
largely analytical means of understanding the results ob- 
tained by these authors. 

In order to describe the zero-temperature phase transi- 
tion from the superfluid to the Mott-insulating phase an- 
alytically, we need to make some appropriate mean-field 
approximation to the hamiltonian in Eq. (|l]). A more or 
less standard approach would be to use the Bogoliubov 
approximation. In Sec. II we show that this approxima- 
tion does not predict the expected phase transition and 
we explain the absence of the phase transition. In Sec. 
Ill we analytically investigate an alternative mean-field 
theory, proposed by Sheshadri et al. and compare 
the analytical results with exact numerical results. In 
Sec. IV we discuss the properties of the Mott insulating 
phases by calculating the quasiparticle and quasihole dis- 
persions and finally in Sec. V we relate the parameters t 
and U to experimental parameters such as laser intensity 
and wavelength. 



II. BOGOLIUBOV APPROXIMATION 

We first transform the hamiltonian to momentum 
space by introducing creation and annihilation operators 
ajj. and ak respectively, such that 



(2) 



t ik-i 
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where A''s is the number of lattice site and is the co- 
ordinate of site i. The wavevector k runs only over the 
first Brillouin zone. For mathematical convenience we 
take only a finite volume V, so that the momenta Kk are 
discretized, which allows us to write sums instead of in- 
tegrals in Eq. (^. Later we will take the continuum limit 
V ^oo. Using the fact that e-'^'^"''') '"' = 7V^(5k,k', it 
is easily shown that the prefactor \/\/Wl ensures that the 
total number of particles obeys N = c\ci = '^^^ aj^Ok- 
If we limit our description to cubic lattices with lattice 
distance a and substitute Eq. into the hamiltonian, 
we find 

= ^ (-ek - aj^flk 

k 

+ ^^mmz!"k4'ak"ak"'4+k',k"+k"', (3) 

k k' k" k'" 

where we defined ek — 2t Yl'j=i cos(fcja), with d the num- 
ber of dimensions. For a Bose condensed gas the average 
number of condensate atoms Nq is a number much larger 
then one, which means that A^o ~ i'^o'^o) ~ (aoio) and 
we are allowed to take iVo = (( 



-.U 



ao • 



Since 



and 



(ao) are complex conjugates, we con- 
clude that (og) ~ (ao) = \/Nq, where we have chosen 
these expectation values to be real. The Bogoliubov ap- 
proach consists of replacing the creation and annihilation 
operators by their average \/Nq plus a fluctuation 



(4) 



and minimizing the energy of the gas with respect to the 
number of condensate atoms A^o- At the minimum, the 
part of the hamiltonian that is linear in the fluctuations 
must therefore be zero. Performing the above subsitution 
and selecting the linear terms yields 



H 



(1) - 



(5) 



where the superscript denotes the order in the fluctu- 
ations. Since H^^^ must be zero for all Aq and oq we 
conclude that in the lowest order approximation 



fi = Uhq — zt, 



(6) 



in terms of the condensate density no = Nq/Ns and the 
number of nearest neighbours z — 2d. This expression 
can be easily understood since the chemical potential 
is the energy needed to add one particle to the system. 
Adding one particle results in an energy increase due to 
the interaction with the no particles already at each site, 
and an energy decrease due to the possible hopping to 
one of z nearest-neighbour sites. 



Next we determine the effective hamiltonian I-f^^ , 
which contains only the parts of zeroth and second or- 
der in the fluctuations. The zeroth-order term is found 
by subtituting all creation and annihilation operators by 
V^/Vq. To find the quadratic term, we substitute in the 
interaction term two creation or annihilation operators at 
a time by \/No and write down all possible combinations. 
Performing the summation over one of the remaining mo- 
menta yields finally 

iJ^ff = (^-2z - M + ^Uno^ No + J2 (-^k - m) flkflk 



k 



(^Okfl-k + 4aj^ak 



"-kftk 



(7) 



We can simplify this expression somewhat by using the 
commutation relation [ak,a]{j.] = 1. If we also substitute 
Eq. (O) and write e^^ — 2z — ek, we find 



(8) 



ek + c/no U no 
U no ek + C/ no 



flk 



where the extra zeroth-order terms are generated by the 
commutation of aj^. and a^. 

The effective hamiltonian is diagonalized by a Bogoli- 
ubov transformation. This implies that we define new 
creation and annihilation operators 6]^ and 5k for which 
the effective hamiltonian is diagonal, by means of 



5k 



Uk Wk 



flk 



B 



flk 



(9) 



To ensure that the operators b^. and 6k still obey the 
standard commutation relations for bosonic creation and 
annihilation operators, we have to demand that the co- 
efficients of matrix B obey 



"k 



«k 



1. 



(10) 



If we now substitute Eq. (||) into Eq. (g) and demand 
that the result reduces to the diagonal hamiltonian 



H'"" = -^UnoNo + IY1 [^^k - (ek + Uno)] 



2 

E 

k 



?iwk5t6k 



(11) 



we find that Uk and f k must be solutions of the following 
two equations: 



((wk)^ + {v].f) Uno - 2ukWk(ek + Uno) = 0, 
(|ukP + |t^k|)^ (ek + Uno) ~ (ukWk + ui^vl)Uno = huj], 



(12) 



Using the normalization in Eq. (p^), we can easily find 
the solution 



2 



(13) 



2 V huJk 



To also obtain the condensate density no, which until 
now has been arbitrary, we now need to calculate the to- 
tal density n as given by our effective hamiltonian. The 
total density is thus given by 



akak)Hrft, 



(14) 



where the brackets ()h<='' denote the expectation value 
as calculated with the effective hamiltonian For a 

Bose condensed gas, this density consists of two parts: 
the density associated with the macroscopic occupation 
of the one-particle ground-state, i.e., the condensate, and 
the density due to the occupation of the higher lying 
one-particle states. In this case, the condensate den- 
sity equals the parameter and the density of the non- 
condensate part is determined by an average over the 
quadratic fluctuations, which will be a function of tiq. 
Calculating the average over the quadratic fluctuations 
by means of Eq. (ph yields first of all 



n = no 



(15) 



k^^O 



If we then use Eq. ( p^ ) and substitute the Bose distribu- 
tion evaluated at hLO\^ for (&j(.&k)_f/t!tf we find that 



n = no 
1 

N. 



(16) 



k^^O 



Ck + Unp 1 



In the zero-tcmperaturc limit, /3 ^ oo, the first term 
in the summant is zero. Taking the continuum limit by 
using V dk/(27r)'^, changing from momenta 

k to q = 27rk/a, and realizing that Ng = V/a'^, we arrive 
at the expression 



no 



1 



1/2 



dq 



-1/2 



(17) 



with 



Eq = 2t^^^Jl - cos(27rgj)] and fiu)^^ = (Cq -I- 
2C/noeq)i/2_ ^e can now obtain the condensate density 
by solving Eq. (^^ for no for a fixed value of n. We ex- 
pect that at integer n, for a fixed value of U /t there will 
be no superfluid solution and this will mark the phase 
transition to the insulating phase as predicted by 



A. Numerical results 

In Fig. 0(a) we plotted the result of this calculation for 
a two dimensional lattice. We see from this figure, that 
there is only a marginal difference between the case that 
n = 0.5 and n = 1.0. In Fig. |](b) we plotted the result 
for a three dimensional lattice. In this case the difference 
between half filling and integer filling is somewhat larger, 
but there is clearly no critical value of U / 1 for which the 
condensate density goes to zero. 

These results lead to the suspicion that the phase tran- 
sition to the insulating phase is not present in this ap- 
proximation. To verify this, we investigate the limit of 
U/t ^ oo \n some detail. 



B. Asymptotic behavior 

When U/t ^ oo we intuitively expect the system to 
become an insulator, because it effectively means that 
the hopping parameter goes to zero. We therefore ex- 
pect that there are no superfluid solutions as U/t — > oo. 
We can see that in this limit the integrand in the right 
hand side of Eq. ( p^ behaves as (C/no/2eq)^/^. One can 
also prove that eq < 47r'^|qpi. This means that 



1/2 



1/2 



dq 



Unp 



> 



+ 2Ue^np 



1 jUnp fi/^ dq 



The integral at the right hand side of Eq. (|18|) can be done 
analytically in two dimensions and numerically in three 
dimensions. When we call the result of the integration 
in d dimensions Li, we see that Eq. (0), for U/t — > oo, 
reduces to 



no 



1 

47r 



Unp _ 1 
2t 2' 



(19) 



where h = 2.22322 and h = 2.38008. This is a quadratic 
equation in ^Jnp which always yields a positive solution 
for no given by 

2 



no 




(20) 



We can correct for the error we made in Eq. ( p^ ) by us- 
ing a higher value for Id, but while this may change the 
value of no, it will still yield a positive solution. We see 
from Eq. ( |20| ) that lim^z/t^oono = as expected, so we 
can conclude that the Bogoliubov approximation as de- 
scribed above does not predict the phase transition to 
the Mott insulator phase in two and three dimensions. 
The reason for this is that the Bogoliubov approach only 
approximately treats the interactions. As a result, the 
Bogoliubov approach cannot describe large depletions of 
the condensate. 

We also see from Eq. (|lj) that in one dimension Ii di- 
verges. Substituting this in Eq. ([T^), we see that there are 



3 



no Bose-condensed solutions ,i.e., solutions with uq 0, 
in one dimension. This is in accordance with the Mermin- 
Wagner-Hohenberg theorem ||To|-[l^ . 

As the Bogoliubov approximation fails to predict the 
phase transition to the Mott insulator phase, we now 
consider a different mean-field theory that treats the in- 
teractions exactly and approximates the kinetic energy 
of the atoms in the optical lattice. 



III. DECOUPLING APPROXIMATION 

To arrive at a mean-field approach, that is capable of 
describing the Mott insulating phase, we start again from 
Eq. (|l|). Analogous to the Bogoliubov approach, we intro- 
duce the superfluid order parameter ip — ^Jnl — (cj) — 
{ci), where rii is the expectation value of the number of 
particles on site i. Note that we take the expectation 
values to be real, as before. We now, however construct 
a consistent mean-field theory by substituting 

c\cj ^ {c^i)c, + c\{c,) - {c\){c,) 

= ^ (c\ + ci) - ^\ (21) 

into Eq. (||). Performing the substitution yields 



i 



(22) 



where z = 2d is again the number of nearest-neighbour 
sites and Ng is the total number of lattice sites, as before. 
This hamiltonian is diagonal with respect to the site in- 
dex i, so we can use an effective onsite hamiltonian. If we 
introduce U — U I zt, fl — ii/zt and the number operator 
■hi = clci, we find 



Hf ^^Uh, (n, 



1) - fxrii 



cj + c, 



(23) 



which is valid on each site i. We will therefore drop the 
subscript i in the following. Note that we scaled all the 
energies by a factor 1/zt, making this hamiltonian a di- 
mensionless operator. 

After writing Eq. (^ ) in matrix form with respect to 
an occupation number basis, we can solve the problem 
numerically by explicitely diagonalizing the part of the 
matrix with occupation number below a certain maxi- 
mum value 1^. Later we also follow this procedure, but 
we first determine the phase diagram analytically using 
second-order perturbation theory. 



A. Second-order perturbation theory 

When we write H""^ = H^"^ + -ipV, with 



= Ic/n (n-1)- p,h + 
F = - (c^ + c) , 



(24) 



we see that in an occupation number basis the odd powers 
of the expansion of the energy in tp will always be zero. 
If we denote the unperturbed energy of the state with 
exactly n particles by eIi ^ , we find that the unperturbed 
ground-state energy is given by 

4°) = to)|n = 0,l,2,...} . . 

^. J mm 

Comparing En'' and E^l^^-^ yields 

(0) ^ r if M < 0, 

^ \^U9{g-l)-fl9 ifU{g-l)<fl<U9. 

(25) 

Next, we calculate the second-order correction to the 
energy with the well-known expression 



\{9\V\n) 



^ p(o) 



E, 



(0) ■ 



(26) 



where |n) denotes the unperturbed wave function with 
n particles, of which the state with n ~ g particles is 
the ground state. Since the interaction V couples only 
to states with one more or one less atom than in the 
ground-state, we find 



^ i7(5-l)-/2 fl-Ug- 



(27) 



If we now follow the usual Landau procedure for second- 
order phase transitions by writing the ground-state en- 
ergy as an expansion in 

EM = ao(5, p) + a2(g, Ji)^^ -f O(V^), (28) 

and minimize it as a function of the superffuid order pa- 
rameter ^, we find that ip = Q when a2{g,U,fL) > 
and that ip 7^ ^ when 02(3, U, p.) < 0. This means that 
0.2(57 U,p,) = signifies the boundary between the super- 
fluid and the insulator phases. Solving 



0-2(9, U,fi) 



g 



Uig~l)-fi p-Ug 



1 = 0, 



yields 



/i± = 2 (C/(2.g - 1) - 1) ± -^C72 - 2C7(25 + 1) + 1, 

(29) 

where the subscript ± denotes the upper and lower halves 
of the Mott insulating regions of phase space. Note that 
this result is exact within our mean-field theory. Fig. 



4 



shows a plot of Eq. ( p9| ) for g = 1,2, 3. By equating 
and fL- we can find the point of smallest U for each lobe. 
Denoting this critical value of U by Uc we have 



Uc = 2g+l + v/(2.g+ 1)2 - 1. (30) 
5.83 for the g = 1 insulator, a value 



which yields C/, 
also found by [f 



B. Fourth-order perturbation theory 

To find out more about the phase transition, we now 
carry out fourth-order perturbation theory to find the 
rate with which the particle density increases as a func- 
tion of /2. In the Appendix, we present a way to calculate 
the higher-order terms in the perturbation series. Using 
this procedure we can write the ground-state energy as 



Eg{i(;) = ao{g,U,fj.) + a2{g ,U , 



ai{g,U,n)^'^ 



(31) 



with 



9(9 - 1) 



{U{g-l)-fl) ([7(25 - 3) - 2/2) 

(g+l)(.g + 2) 
{fi-Ugf{2fl~U{2g + l)) 
f 9 , 5+1 



U{g-l)~ii fi-Ug 



1 



■ (32) 

\{Uig~l)-fiy {fl-UgYJ 

In Figs, ^(a) and (b) we show plots of Eq. ( ^ together 
with the result of an exact numerical diagonalization of 
the effective hamiltonian. As can be seen, the overlap is 
very good near the boundary given by Eq. (p9|). In Fig. 
^ it can be seen that the numerical result exhibits a cusp 
when U = fl, which is not predicted by Eq. This 
is due to the fact that in this particular case we need to 
use first-order degenerate perturbation theory, because 
at fi — nU the states with n — 1 and n particles per site 
form a doubly degenerate ground-state. The resulting ex- 
pression for the ground state energy is now nonanalytic 
and given by 



1 



f_L—nU 



Un{n+l) + ij^ (33) 



which is the solid line in Fig. |^. Note that the occurence 
of a cusp is analogous to the well-known Jahn- Teller ef- 
fect in solid-state physics 

We now continue by calculating the average number of 
particles per site in the grand-canonical ensemble by 



d{H 



cff\ 



9 - ^ 



dfj, 

d_ / a2{g,U,jif 
dp. \ Aa4{g,U,fl) 



(34) 



where ■0min = [— ^2(3, C7, M)/2a4(g, C7, /l)]^/^ jg the mini- 
mum of Eq. (|3l|). Making use of the previous results, we 
can now plot the density as a function of p, for a fixed 
value of U. Between the edges p±, the density will re- 
main constant because V'min — and the second term in 
the right-hand side of Eq. (|3^ ) does not contribute. Out- 
side that region, the density will start to change with a 
nonzero slope. In Fig. ^ this is plotted for U = 11. The 
solid line shows the result of the calculation described 
above and the dash-dotted line is a numerical result ob- 
tained by exact diagonalization. As can be seen, the ana- 
lytical results are in good agreement with the numerical 
calculation. We can now also qualitatively understand 
the solution found numerically by Jaksch et al. ^ for an 
optical lattice in an external harmonic trap. In a first ap- 
proximation, we can describe the effect of a slowly vary- 
ing trapping potential by substituting p, in Eq. (0)by 
p' = p + V(r), where V{r) is the external trapping po- 
tential. Combining this with Fig. || yields the density 
profile found in P|. 



IV. DISPERSION RELATIONS 



An important property of the Mott insulating phase is 
that the fluctuation in the average number of particles 
per site goes to zero at zero-temperature. Since these 
fluctuations can be described as quasiparticle and quasi- 
hole excitations, we will study these now. We calculate 
the quasiparticle and quasihole dispersions using a func- 
tional integral formalism. We start by deriving an ex- 
pression for the effective action. Readers unfamiliar with 
functional integrals may want to skip to subsection IV B, 
where we discuss the results of the calculation. 



A. The effective action 



We define complex functions a*(T) and ai{T), respec- 
tively, and write the grand-canonical partition function 
as 

Z = Tre"^^= I Va*Vaexp{~S[a*,a]/fi}, (35) 



where the action S[a* , a] is given by 



S[a* , a] — dr 
Jo 



(36) 
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with j3 = l/fcsT, ks Boltzmann's constant and T the and integrate over the original fields a* and a^. If we 
temperature. To decouple the hopping term, we per- 
form a Hubbard- Stratonovich transformation by adding plicitely that 
a complete square to the action, which then becomes 



S[a* , a, -0* , -0] — S[a* ,a 
f.n(3 



(37) 



"'0 



Here ip* and ip are the order parameter fields. To obtain 
an effective action as a function of these fields, we rewrite 
Eq. @ 



as 



denote by 5"*^°^ [a*, a] the action for — 0, we have ex- 

exp {-S'''[r,^]/ri) = exp ^-i jj ^^Y. ^'W-'^V'. j 
X j 2?a*2?acxp|-S'(°)[a*,a]/?i| 



X exp 



S[a*,a,tlj*,tp] = 



dr 



(38) 

We can now calculate 5'°'^ perturbatively by Taylor 
expanding the exponent in the integrant of Eq. ( |39| ) and 
evaluating the various correlation functions of the field 
, theory given by S^^\ This yields for the quadratic part 
of the effective action 



\ \ y / / g(o) u 



1 

m 



"'0 



iji'j' 



/ dT^Ujil^li^j. 

■JO 



(40) 



If we perform the multiplication in the first term in the 
right-hand side and use the information we have about 
the correlations in the unperturbed system, i.e., 

«aj*)5(o) = (aiaj>s(o) = 0' 

(a*aj)s(o) = 3(0) = (oiOsw ^hJ^ (41) 

we obtain in first instance 
phu ( 

5(2) [r,^]= drlY, U,rdr)Ur) (42) 
~ % j X! ^»J*«'i'V'j(T) (ai(r)a*,(T'))^(o) i/'/(^') 



where we have now shown the r dependence of the fields 
explicitely for clarity reasons. Because we will only con- 
sider nearest-neighbour hopping, we write 

_ _ J i for nearest neighbours 
t^j - i'H - \ otherwise. ^ > 



First we treat the part of Eq. (|4^) that is linear in tij . 
We have 

E (^)^J- (^) = E (^)^»±{i} (^)' (44) 



where ±{1} denotes all possible jumps to nearest neigh- 
bours. In the case of one dimension this would simply be 
±1. If we call the lattice spacing a and introduce carte- 
sian momentum components ki with i = 1, • • • , d, where 
d is again the number of dimensions, we find 

d 

E*ijV'r(T)V'j(T) = E2^^k(T)0k(^)E*=°^(^J")■ 
^j k j=l 

(45) 

Next we calculate the part that is quadratic in tij. We 
can treat this part by looking at double jumps. The 
expectation value of {aiali) gi^o) is proportional 5ni and 
independent of the site i, according to Eq. (|4l|). This 
means that we find, with similar notation as before, 

E *y^i'j''^i(^) (a«(^X'(^')>sW) V'j'(t') 

ji'j' 

= ("^(^x (t'))s(o) E^y V'^j(^)'^j'(^') 

jj' 

= e (a,(r)<(r'))<,<o,E{^^;(^)V-'.(^') 
j 

+ v^;(t)V',±{2}(t') + V';(T)v,±{y2}(^')} , (46) 
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with z again the number of nearest neighbours. The first 
term in the summant is a jump in each direction, fol- 
lowed by a jump back. The second term indicates two 
jumps in the same direction and the third term is a jump 
in each direction followed by a jump in a perpendicular 
direction. Note that the third term is absent in one di- 
mension. It can be shown that the complete double jump 
term reduces to 

ji'j' 

= {adT)aUr'))sio, E ^kWV'k(r')(ek)^ (47) 



where we again used ek = '2tJ2j=i cos{kja). 

To also treat the time dependence of the fields, we in- 
troduce Matsubara frequencies huJn = iT{2n)/hf3 by 



n k 



dr'gexpl- {ihtJn + IJ.- {g - 1)U)t' /h} 



(53) 



Performing the r' integration we then easily obtain 
5(2)[^*,^]=^^|V,k„|2(ek)x (54) 

9 



l-(£-k) 



71 k 

3 + 1 

1 

-ihuJn — M + 9U ihuJri 



1)U 



(48) 



To translate the expectation value of the fields into the 
expectation value of operators, we introduce an (imagi- 
nary) time ordering operator T. As a result 



Note that this result is exact within our mean-field 
theory. It contains all powers of the frequencies and 
momenta and no gradient expansion has been applied. 
This is important because the elementary excitation are 
gapped as we will show in the next section. 

We can find an equation for real energies huj by sub- 
tituting iujn — > and equating the remaining factor to 
zero. This gives 



: 



l-(fk) 



5 + 1 



-tiu; — fJ, + gU huj + ji — [g — 1)U 



(55) 



(ai(T)a*,(T'))5(o) 



T 



S(0) 



(49) 



The time ordering can easily be expressed in Heavyside 
functions as 



T 



«.(-)4(-')])^,„, =6{r-r'){adr)4{r'))^^^^ 

' ~r)ial{T')adr))^^^^. (50) 



Ultimately this yields the result Eq. ( |56| ) given below in 
subsection IV B. 



B. Results 

Now we will explore the results of the calculation pre- 
sented in the previous subsection. The quasiparticle and 
quasihole dispersions are given by 



If we use the unperturbed energies as given by (|2^), we 
thus find 

(ai(r)a*,(r'))5(o) 
= e{T - r')(l + g) cxp {- (4% - 4")) (r - T')/h] 

+ 6{r' - r)gexp { (Ef\ - Ef^) (r - r')/n]. (51) 

Because g minimizes E^'^ we know that 
Efl,-Ef^ = -^^ + gU>0, 



tiuj, 



qp,qh 



1 



± 2V(e"k)'-(45 + 2)C/ek + f/2. 



(56) 



In Fig. ||(a) we show for k = a plot of the above equa- 
tions. The dotted lines indicate the asymptotes of Eq. 
(^), which are given by 

lim hujqp = + gU — (g + l)eo 

f7— »oo 



E. 



(0) 

9+1 



E^°^-{g + l)zt, 



(52) 



lim huj, 

U^oo 



qh 



-fi+{g-l)U + geQ 



4«)-i?fi = ~^-t-(.g-l)C/<0. 



= 4") - <\ + gzt, 



(57) 



Note that we use parameters /i and U instead of p, and 
U, because we have not yet divided out the factor zt. 
Combining the above with Eq. (E2|) wc find 



with - E'g"' and E'g"' - £;^_\ given by Eq. (^). The 
difference between Eq. (|5^) and Eq. (|5^) is caused by the 
fact that Eq. (m) is calculated for t = 0. It can easily 
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be understood that for t 7^ 0, the first-order correction 
is due to the hopping terms c^c^t, where site j is one of 
the nearest neighbours of site i. When we have g parti- 
cles in all lattice sites and we add one particle to site i, 
we have (cjci) = g + 1, so the effective hopping parame- 
ter for a particle is {g + l)t. However, when we remove 
a particle from site i, we have (c|cj) = g, which repre- 
sents a particle hopping to site i from one of its nearest 
neighbours. The effective hopping parameter for a hole 
is therefore only gt. In combination, we see that in the 
limit of [/ ^ 00, Eq. ( |56| ) indeed reduces to a physically 
intuitive result. 

As shown above, the slopes of the asymptotes differ 
exactly by U, so in the limit of U / zt — > 00 the gap for 
the creation of a quasiparticle-quasihole pair is equal to 
U . We can find a first approximation for the dispersion 
of the density fluctuations by substracting the two solu- 
tions, which yields 



qp 



hujqh 



V'(ek)2-(4g + 2)?7ek + t/2. 



(58) 



In Fig. ^(b) we show again for k = a plot of the above 
equation as a function oi tj — U j zt for g = 1. We can 
see that there is a band gap, which proves that the MI 
phase is indeed an insulator and we also see that the 
band gap disappears as we approach the critical value 
Uc = Uc/zt K, 5.83 that was found earlier. For smaller 
values of U we are in the superfluid phase, which ac- 
cording to the Hugenholtz-Pines theorem is expected to 
always have gapless density fluctuations. 



V. MICROSCOPIC PARAMETERS 

To estimate the experimental feasibility of observing 
the described phase transition, we now relate the param- 
eters t and U to the microscopic parameters. Because we 
have an experimental interest in sodium, we will calcu- 
late these parameters for sodium atoms trapped in a lat- 
tice constructed with four laser beams. To calculate the 
hopping parameter t, we calculate the overlap between 
single particle ground-state wave functions in neighbor- 
ing sites. To calculate the interaction strength U , we use 
the pseudopotential method. 

First we calculate the optical potential experienced by 
the atoms, following the approach of Petsas et al. [ [14[ . 
We describe a J = 1/2^ J = 3/2 transition and choose 
a laser beam conflguration with two pairs of laser beams. 
Each pair lies in a plane and the planes are perpendic- 
ular to each other. All beams have the same angle 
with respect to the intersection of the two planes. We 
choose the quantization axis along the intersection and 
label it as the z-axis. Furthermore, we choose the po- 
larization of the laser beams linear and perpendicular to 
the plane spanned by the pairs of beams. The conflgu- 
ration is illustrated in Fig. 0. It should be noted that it 



is also possible to simply superimpose d standing waves 
to obtain a d dimensional lattice, but this requires sta- 
bilization of the relative phases of the laser beams. We 
define /f, as the sum of the intensity of the beams and if 
k is the magnitude of the k- vector, we define k± — k sin0 
and kj I —k cos d. If we add the electric field components 
and express them in spherical components 



iE. 



Eq — Ez, 



(59) 



we find that the spatial dependence of the intensity of 
the resulting light field is given by 

I±/h = ^ (cos^ {k±x) + cos^ (fc_Ly) 

± 2 cos (fc_Lx) cos (A:_Lj/) cos (2fc//2:)) , 

h/h = 0. (60) 

Note that at the minima of I± the polarization is purely 
cr*. Also note that since the linear component is always 
zero, the two ground-state levels are not coupled. 

Following Nienhuis et al. [|l^ we can now calculated the 
optical potential. Because of the fact that the ground- 
states are not coupled, we can treat them seperately. 
With 5 the detuning, F the rate of spontaneous decay and 
Q.± the Rabi frequencies for the cr* components of the 
light field, we can write the potential for the rrij = ±1/2 
level in the limit of low saturation as: 



V± 



1 



h5 



2 1 + 4((5/r)2 



21174 



12|f]=p|2 



3 F^ 



(61) 



where the factor ^ arises because of the Clebsch-Gordan 
coefficients for J = 1/2 ^ J = 3/2 transitions. 
Now we define a convenient prefactor: 



hSsQ 



2 1 + 4((5/F)2 



h6s 



(62) 



where s = sq/ (l + 4((5/r)^) is the off-resonance satura- 
tion parameter and so — 2|rip/F^ is the on-resonance 
saturation parameter, which is usually written as sq = 
I /Is- The saturation intensity Is is a constant for a given 
transition. If we substitute Eq. (|62[) in Eq. (|6l|), we find 



V± = Vb 



h 3 h 



= (cos^ {k±x) + cos^ {k±y) 
± cos (/cj^x) cos (fc_Ly) cos (2A://z)) 



(63) 



We now write the hamiltonian of a particle in the po- 
tential as: 
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opt 



2m 



(64) 



and solve the time-independent Schrodinger equation 
variationally by assuming an isotropic gaussian wave 
function and minimizing the energy as a function of width 
of the gaussian. If we call (3 the width of the wave func- 
tion, we can write the normalized wave function as 



*(r) 



1 



3/4 



(65) 



We assume we have a spin polarized sample of atoms, so 
we can use either the V+ or the V- potential. 

For simplicity we now calculate the parameters for a 
one dimensional lattice. For the lattice configuration in 
Fig. this gives approximate results, but for a phase sta- 
bilized superposition of three standing waves, the results 
are immediately applicable. In this case, the potential 
reduces to V± = 2H[2 ± cos(2/cz)]/3 + ki^{x'^ + V^)l'^. 
where the transverse potential is caused by the fact that 
the laser beam has a finite width. If we assume the wave- 
function is tightly localized in the center of the local po- 
tential well, we can approximate the potential as an har- 
monic potential V+ = 214 + Kr^/2 with k = — 8Vf,fc^/3, 
where we assume we can adjust the width of the laser 
such that A«_L « K. These approximations yield the well- 
known equations for the width and the level splitting in 
the potential 



/3= — 

TYIK 



1/4 



(66) 



Using the above width (3 we calculate the value of the 
interaction strength U with the pseudopotential method. 
According to this is valid for sodium even if the width 
of the trapping volume is of the order of the scattering 
length. In general the interaction strength between two 
atoms in the same one particle wave function is given by 

U = j dr j dr'**(r)**(r')M„t(r- r')^'(r)^'(r'), (67) 

where 14nt(r — r') is the interaction potential. If we ap- 
proximate the potential as 



^nt(r-r') = 



^(r-r'). 



(68) 



we can write Eq. (|6^) as 
At: a. It? 



U = 



m 
m 

2TiLo ( a 



dr^'*(r)**(r)*(r)*(r) 



dr|*(r) 



(69) 



where is the triplet s-wave scattering length. Accord- 
ing to the value of the scattering length for a spin 



polarized sodium-sodium collision is Os = (85 ± 3)ao. 
Note that the use of a one-band model is justified when 
?7 < or /3 > 2a,,/(27r)i/2 ~ 3. 

onm. 

Next we calculate the value of the hopping parameter 
t. In the tight binding limit t is given by 



i= - / dr**(r) ( ^ + Kt ) *(r + aej), (70) 



where ij is an axis vector along a lattice direction, so that 
when ^'(r) is the ground-state wave function, ^'(r-|~aej) 
is the ground-state wave function of an atom at a neigh- 
boring site. One can show that product of two wave- 
functions at neighbouring sites is a gaussian function 
centered around r -I- aej/2. We can therefore approxi- 
mate the potential around the maximum of the barrier 
by: V± = 2V6/3T KrV2. Substituting this into Eq. (|7C 
ultimately yields 



t = 



1 - 



(71) 



Figs. ||(a) and (b) show plots olU /Ej. and t/ E,. respec- 
tively, with Er ~ h^k^/2m the recoil energy. Both are 
plotted as a function of the trap depth T4rap = — 4Vb/3. 
The values were calculated for a laser wavelength of 600 
nm. The saturation parameter needed to reach these trap 
depths is in the order of 10^, which is not imrealistic ex- 
perimentally. 

Fig. ^ shows also U/zt as a function of the trap depth, 
for a wavelength of 600 nm, in one, two and three dimen- 
sions. Again, the saturation parameter is in the order 
of 10^. As can been seen, the desired critical value is 
reached in all three dimensions. The value of the width /3 
lies between 12% and 8% of the wavelength in the range 
considered in the above plots. This implies that both 
the harmonic approximation and the use of the one-band 
model are justified. 



VI. CONCLUSION 



Due to the absence of the superfluid-insulator phase 
transition in the Bogoliubov approach, we conclude that 
the interaction is the dominant component in this phase 
transition. When the interaction energy is treated ex- 
actly, the theory indeed predicts a phase transition. The 
mean-field theory predicts a phase transition even in one 
dimension, which we expect to survive as a Kosterlitz- 
Thouless transition when fluctuations are incorperated. 
However a definite prove of this requires further study. 

We analytically calculated the phase diagram and the 
particle and hole dispersion relations in the insulator 
phase. A first-order approximation to the dispersion of 
the density fluctuations shows that the system indeed 
goes from a gapped to a gapless phase. A calculation of 
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this dispersion below the critical value for 17/ zt will have 

to be done in order to check the presence of linear dis- = ^ (a|y|n) 
persion that would verify the assumption that the phase n,p,q^a 
with ■!/) is indeed superfluid. The one-band model we 
used to calculate the parameters for sodium gives opti- 
mistic results for future experiments, within the range of 



_^(2) {n\V\a) 



( 

-^a'-T— 

\y\v) WW) w\<^) 



(0) _ p(0)\2 



parameters it allows. {^'a ^ ^ E^^^ ^Ea — Ep -'^ (^si 



(0) p(0) 



APPENDIX: THE PERTURBATION SERIES 



(A5) 



A powerful approach to calculating higher-order terms 
in the perturbation series is derived in [[l8[ . Here we only 
give the result of that derivation. If we denote by \n) the 
unperturbed wave functions and i?!*''' the unperturbed 
energies, we can define an operator 



S. 



-\a){a\ 



if fc = 



si"') 



and one can prove that the n-th order correction on the 
energy e'^^ is given by 



£;(») ^ Tr 



{n-l} 



(A2) 



where {n} = {fco) ^n|fco + •■• + fcn = In the case of 
n = 2, this quickly gives the well-known result 



= Tr 



{1} 

a\VSlV\a), 



E 



{n\V\a) 



— I 



AO) e(o) 



(A3) 



The same can be done for Ei^\ ' and si'^K The first 
two can be shown to involve only terms proportional to 
odd orders of V and with V gc (c^ + c) these are of course 
zero. The fourth-order term is in general given by 



.{3} 

= {a\VSlVSlVSlV\a) - {a\VSlV\a){a\VSlV\a) 
-2{a\V\a){a\VSlVSlV\a) + {a\V\af {a\VSlV\a) . 

(A4) 

If we drop the terms containing expectation values of odd 
powers of V and substitute Eq. (Al), we find 



(Al) which we have used to derive Eq. (32) 



Ei''^ = {a\VSlVSlVSlV\a) - {a\VSlV\a){a\VSlV\a) 
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FIG. 1. The condensate fraction no/n (a) in a 
two-dimensional optical lattice and (b) a three-dimensional ' ' ' ''° 

^ ' -0.50 -0.25 0.00 0.25 0.50 

optical lattice, both as a function of U jt for n = 0.5 (dashed if 
line) and n = 1.0 (dotted line). 




FIG. 2. Phase diagram of the Bose-Hubbard Hamiltonian 
as obtained from second-order perturbation theory (solid 
lines). The dotted lines indicate the zeroth-order phase di- 
agram. Later on, Fig. ^ is taken along the dashed line in this 
figure. 



FIG. 3. Ground-state energy as a function of i/; for (a) 
f7 = 11 and = 8.9 and for (b) f7 = 11 and p = 7.8. The solid 
line represents fourth-order perturbation theory whereas the 
dotted line represents a numerical diagonalization of the ef- 
fective hamiltonian. The dashed line is the difference between 
the two (scale on the right). 




(a) 




3.910 [ 
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0.00 



] -0.020 

0.50 



FIG. 4. Ground-state energy as a function of il) for 
f7 = /i = 11, as obtained from first-order perturbation theory 
(solid line) and from numerical diagonalization of the effective 
hamiltonian (dotted line). 
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FIG. 5. The density as a function of chemical potential 
= ^/zt for an interaction strength of 17 = U / zt — 11, i.e.. 



along the dashed line in Fig. ^. 
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U/zt 



FIG. 6. The quasiparticle and the quasihole energy for 
k = in the g = 1 insulator lobe. The dotted lines are 
the asymptotes of the curves. The inset shows the result- 
ing first-order approximation to the dispersion of the density 
fluctuations. 




FIG. 7. Laser beam configuration for a three-dimensional 
optical lattice. 
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FIG. 8. Plot of (a) U and (b) t as a function of the trap 
depth. All quantities are in units of the recoil energy Er. 
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FIG. 9. The dimensionless parameter U/zt plotted as a 
function of the trap depth for 1 (solid line), 2 (dotted line) 
and 3 dimensions (dashed line). The dash-dotted line is the 
critical value U = 5.83. 
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